﻿#region Copyright information
// 
// Copyright © 2005-2013 Yongkee Cho. All rights reserved.
// 
// This code is a part of the Biological Object Library and governed under the terms of the
// GNU Lesser General  Public License (LGPL) version 2.1 which accompanies this distribution.
// For more information on the LGPL, please visit http://bol.codeplex.com/license.
// 
// - Filename: GeometricDistribution.cs
// - Author: Yongkee Cho
// - Email: yongkeecho@gmail.com
// - Date Created: 2012-09-06 11:39 AM
// - Last Modified: 2013-01-25 3:59 PM
// 
#endregion
using System;
using System.Collections.Generic;
using System.Linq;

namespace BOL.Maths.Distributions
{
    /// <summary>
    /// Represents a geometric distribution with probability parameter p.
    /// <see>Mathematica: ref/GeometricDistribution</see>
    /// <see>http://en.wikipedia.org/wiki/Geometric_distribution</see>
    /// </summary>
    public sealed class GeometricDistribution : IUnivariateDistribution<int>, IEquatable<GeometricDistribution>
    {
        #region private variables

        private double _p; // the probability of success
        private double _q; // the probability of failure

        #endregion

        #region Public properties

        /// <summary>Gets or sets the probability of success of a geometric distribution.</summary>
        public double P
        {
            get { return _p; }
            set
            {
                if (value <= 0 || value > 1)
                    throw new ArgumentOutOfRangeException("value", "p must be in interval (0, 1].");

                _p = value;
                _q = 1 - _p;
            }
        }

        /// <summary>Gets or sets the probability of failure of a geometric distribution.</summary>
        public double Q
        {
            get { return _q; }
            set
            {
                if (value <= 0 || value > 1)
                    throw new ArgumentOutOfRangeException("value", "q must be in interval (0, 1].");

                _q = value;
                _p = 1 - _q;
            }
        }

        /// <summary>Gets the range of a geometric distribution.</summary>
        public IRange<int> Domain { get { return new Range<int>(1, int.MaxValue); } }

        /// <summary>Gets the mean of a geometric distribution.</summary>
        public double Mean { get { return 1 / _p; } }

        /// <summary>Gets the median of a geometric distribution.</summary>
        public double Median { get { return -Math.Ceiling(Math.Log(2) / Math.Log(_q)); } }

        /// <summary>Gets the mode of a geometric distribution.</summary>
        public double Mode { get { return 1; } }

        /// <summary>Gets the variance of a geometric distribution.</summary>
        public double Variance { get { return _q / (_p * _p); } }

        /// <summary>Gets the skewness of a geometric distribution.</summary>
        public double Skewness { get { return (2 - _p) / Math.Sqrt(_q); } }

        /// <summary>Gets the kurtosis of a geometric distribution.</summary>
        public double Kurtosis { get { return 6 + _p * _p / _q; } }

        /// <summary>Gets the entropy of a geometric distribution.</summary>
        public double Entropy { get { return (-_q * Math.Log(_q, 2) - _p * Math.Log(_p, 2)) / _p; } }

        #endregion

        #region Constructor

        /// <summary>
        /// Instantiates a geometric distribution.
        /// </summary>
        /// <param name="p">probability of success</param>
        public GeometricDistribution(double p)
        {
            if (p <= 0 || p > 1)
                throw new ArgumentOutOfRangeException("p", "p must be in interval (0, 1].");

            _p = p;
            _q = 1 - p;
        }

        #endregion

        #region ICloneable implementation

        public GeometricDistribution Clone()
        {
            return new GeometricDistribution(_p);
        }

        object ICloneable.Clone()
        {
            return Clone();
        }
        
        #endregion

        #region public methods

        /// <summary>
        /// Returns the geometric density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>geometric density</returns>
        public double Pdf(int value)
        {
            if (value < 0)
                throw new ArgumentOutOfRangeException("value", "value must be a positive integer.");

            return Math.Pow(_q, value) * _p;
        }

        /// <summary>
        /// Returns the cummulative geometric density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>cummulative geometric density</returns>
        public double Cdf(int value)
        {
            if (value < 1)
                throw new ArgumentOutOfRangeException("value", "value must be a positive integer.");

            return 1 - Math.Pow(_q, value + 1);
        }

        /// <summary>
        /// Returns the value of the geometric distribution for given probability.
        /// </summary>
        /// <param name="p">probability</param>
        /// <returns></returns>
        public int Quantile(double p)
        {
            return (int)(Math.Log(1 - p) / Math.Log(_q));
        }

        /// <summary>
        /// Returns a random value from the geometric distribution.
        /// <see>http://people.csail.mit.edu/milch/blog/index.html</see>
        /// <see>Reference: Luc Devroye, Non-Uniform Random Variate Generation, p. 87 (available at http://cg.scs.carleton.ca/~luc/rnbookindex.html)</see>
        /// </summary>
        /// <returns>random value</returns>
        public int Sample(Random r)
        {
            return (int)(Math.Log(r.NextDouble()) / Math.Log(_q));
        }

        /// <summary>
        /// Estimates parameters of the geometric distribution using maximum likelihood.
        /// <see>http://en.wikipedia.org/wiki/Geometric_distribution#_parameter_estimation</see>
        /// </summary>
        public void MaximumLikelihoodEstimate(IEnumerable<int> source)
        {
            P = 1 / (1 + source.Average());
        }

        public void BayesianEstimate(IEnumerable<int> source, params dynamic[] priorParameters)
        {
            throw new NotImplementedException();
        }

        #endregion

        #region IEquatable<GeometricDistribution> implementation

        public bool Equals(GeometricDistribution other)
        {
            return _p.Equals(other._p);
        }

        #endregion

        #region Object overriden

        public override int GetHashCode()
        {
            return P.GetHashCode();
        }

        public override bool Equals(object other)
        {
            if (other == null)
                throw new ArgumentNullException("other");

            if (!(other is GeometricDistribution))
                throw new InvalidCastException("The 'other' argument is not a GeometricDistribution object.");
            
            return Equals(other as GeometricDistribution);
        }

        public override string ToString()
        {
            return String.Format("Geometric(p = {0:G4})", _p);
        }

        #endregion
    }
}
